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HIGHLIGHTS 


►  A  distributed  thermal  model  for  a  lithium-ion  electrode  plate  pair  is  developed. 

►  This  model  has  multiple  dimensions  and  multiple  length  scales. 

►  The  model  was  developed  by  coupling  the  heat  equation  with  a  P2D  electrochemical  model. 

►  Reduced  order  model  was  also  produced  to  significantly  reduce  the  simulation  time. 
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This  paper  presents  a  distributed  thermal  model  for  a  lithium-ion  electrode  plate  pair  used  to  predict  the 
distributed  electrical  and  thermal  behavior  of  the  electrode  pair  including  tabs.  Our  model  was  devel¬ 
oped  by  coupling  the  heat  equation  with  a  pseudo  two  dimensional  (P2D)  physics-based  electrochemical 
model.  The  local  heat  generation  rate  is  predicted  by  the  P2D  model  at  every  node  point  in  the  2D 
electrode  pair.  To  reduce  significantly  the  computation  load  of  the  model,  a  linear  approximation  method 
is  introduced  to  decouple  the  electrochemical  model  from  the  heat  equation  with  a  very  slight  loss  in 
accuracy. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Currently,  large-format  Li-ion  cells  are  widely  used  in  the  hybrid 
electrical  vehicles  (HEV),  and  these  cells  are  often  operated  at  very 
harsh  electrical  and  thermal  conditions  such  as  the  high  and  fast¬ 
changing  current  rates  or  extreme  ambient  temperatures.  Under 
such  conditions,  the  thermal  and  electrical  behavior  will  be  quite 
non-uniform  through  the  cell  volume  [1-4].  Multi-scale  and  multi¬ 
dimensional  (MSMD)  modeling  approaches  have  been  proposed  to 
simulate  the  distributed  thermal,  electrical,  and  chemical  behavior 
of  large  format  Li-ion  cells  [5-7].  The  length  scales  of  computa¬ 
tional  sub-domains  in  an  MSMD  model  have  typically  three 
different  levels:  the  microscopic,  mesoscopic,  and  macroscopic 
[7,8].  Modeling  in  the  microscopic  scale  generally  focuses  on  the 
molecular-level  quantum  behavior,  modeling  in  the  mesoscopic 
scale  specifically  deals  with  the  electrochemical  processes  and 
transport  phenomena  in  the  electrode  coatings  and  the  active 
material  particles,  and  modeling  in  the  macroscopic  scale  mainly 
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describes  the  non-uniformity  of  the  electric  potential  along  the 
current  collectors  in  electrode  composites  and  the  temperature 
throughout  the  cell  volume  [7,8].  In  the  fields  of  engineering  and 
industry,  models  with  mesoscopic-to-macroscopic  length  scales 
are  more  preferable  for  application;  and  among  these  models,  the 
electrochemical-thermal  (ECT)  coupled  models  [1]  are  commonly 
used  for  the  thermal  management  design  and  thermal  runaway 
analysis. 

Unfortunately,  the  physics-based  MSMD  models  require 
a  significant  amount  of  computation  time.  The  widely  accepted 
physics-based  model  for  Li-ion  cells,  Newman’s  pseudo-2D  (P2D) 
porous  electrode  model  [9-11],  involves  a  large  non-linear  DAE 
system  in  the  mesoscopic  scale.  Although  various  model  order 
reduction  techniques  have  been  suggested  to  simplify  the  New¬ 
man’s  P2D  model  and  reduce  the  computation  time  for  a  MSMD 
model  [12-17],  the  reduced  order  models  (ROM)  still  have  disad¬ 
vantages  in  terms  of  the  accuracy  for  high  current  rates  simulation. 
For  example,  the  SVM  model,  which  was  derived  by  Smith  et  al. 
[12,13]  and  employed  in  the  MSMD  model  by  G-H  Kim  et  al.  [7] 
show  poor  agreement  with  the  full-order  model  in  the  voltage 
response  when  the  current  rate  exceeds  5C.  Therefore,  the  goal  of 
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good  time-efficiency  in  simulation  as  well  as  sufficient  accuracy  in 
predictions  is  a  great  challenge  for  researchers  in  the  MSMD 
modeling  for  the  Li-ion  cells. 

In  this  work,  we  present  our  full  physics-based  (P2D)  MSMD 
model  and  an  associated  ROM  for  a  Li-ion  electrode  pair  with 
planar  electrode.  We  derived  our  ROM  by  segregating  the  meso¬ 
scopic  and  the  macroscopic  sub-domains,  that  is,  we  decoupled  the 
Newman’s  P2D  model  from  the  charge  balance  in  the  current 
collectors  and  the  energy  balance  in  the  cell  volume.  Significant 
improvements  in  the  simulation  speed  were  achieved  through  our 
model  order  reduction  approach.  To  ensure  that  our  ROM  also 
maintains  sufficient  accuracy  in  predictions,  extensive  model-to- 
model  validations  were  made  between  the  ROM  and  the  full- 
order  full-distribution  model. 

2.  Mathematical  model 

In  this  part  of  work,  the  model  developments  and  simulation 
case  studies  are  based  on  one  single  electrode  plate  pair;  and  in  the 
future,  the  model  will  be  extended  to  cell  stacks  that  include 
multiple  electrode  pairs. 


Table  1 

The  dependent  variables  in  different  sub-domains. 


Length  scale 

Sub-domain 

Dependent  variable 

Symbol 

Independent 

name 

variables 

Macroscopic 

Cell  and  tabs 

Temperature  (K) 

T 

tXY 

Current  collector 
potential  (V) 

$/(/  =  P,n) 

tXY 

Mesoscopic 

Porous 

Li  concentration  in 

ce 

t,xXY 

electrode 

electrolyte  (mol  rrr3) 
Potential  in 

02 

t,xX,Y 

solution  phase  (V) 
Potential  in  solid 

01 

t,xX,Y 

phase  (V) 

Surface  reaction 
rate  (mol  m-2  s-1) 

JjU  =  p,n) 

t,xX,Y 

Particle 

Li  concentration  in 
solid  phase  (mol  m-3) 

cSj(j  =  P,n) 

t,r,xX,Y 

presented  in  Fig.  2.  The  governing  equations  for  the  charge  balance 
in  the  current  collectors  of  the  cell  sub-domain  are  as  follow: 


a  -—l  +  a  + 

CCJ9X2  +  CCJ  dY2  +  6j 


(1) 


2.1.  The  modeling  scales  and  sub-domains 

The  planar  electrode  plate  pair  of  the  Li-ion  cell  can  be  divided 
into  several  computational  sub-domains  in  terms  of  length  scales 
(see  Fig.  1 ).  As  shown  in  Fig.  1,  the  thickness  of  electrode  pair  (the  x 
dimension)  and  the  radius  of  particles  (the  r  dimension)  are  much 
smaller  than  the  length  and  width  of  the  electrode  pair  (the  X  and  Y 
dimensions).  Therefore,  in  our  model,  the  macroscopic  sub- 
domains  include  the  electroactive  part  of  the  electrode  plate  pair 
(the  cell  sub-domain)  and  the  two  electrode  tabs;  the  mesoscopic 
sub-domains  include  the  porous  electrode  coatings  and  separator 
(the  porous  electrode  sub-domain),  and  the  active  material  parti¬ 
cles  (the  particle-sub-domain).  According  to  a  dimensional  anal¬ 
ysis,  the  temperature  and  the  electrical  potential  of  each  current 
collector  can  be  considered  to  be  distributed  only  in  the  two 
macroscopic  (X  and  Y)  dimensions;  the  mass  and  charge  transport 
in  the  electrolyte  and  the  solid  phase  in  the  porous  electrode  sub- 
domain  follow  the  mesoscopic  (x)  dimension;  the  mass  balance  in 
the  particle  sub-domain  is  described  by  another  mesoscopic 
dimension,  the  r  dimension.  Table  1  summarizes  the  dependent 
variables  to  be  solved  in  the  different  sub-domains. 

2.2.  The  models  in  the  macroscopic  sub-domains 


where  <f>j  is  the  potential  in  the  current  collector  of  electrode  j,  accj 
is  the  electrical  conductivity  of  current  collector,  8j  is  the  thickness 
of  current  collector,  and  i^j  is  the  normal  inward  current  density 
through  the  coating/current-collector  interfaces.  In  the  tab  sub- 
domains,  there  is  no  transverse  current,  and  the  charge  balance 
equations  are: 


a  2<f>j  a2^- 

<Tccj~dXT  +  <Tccj~dYT 


o  (j  =  P,n) 


(2) 


The  boundary  conditions  (see  Fig.  2)  at  the  tops  of  electrode  tabs 
are: 


n'(  *ccjVfc/)  —  lapp  j 


(3) 


where  n  is  the  unit  normal  vector  pointing  out  of  the  boundary,  and 
iappj  is  the  applied  outward  current  density  based  on  the  cross 
section  at  the  top  of  each  electrode  tab. 


2.2.2.  The  energy  balance 

The  equation  for  the  energy  balance  in  the  cell  and  tab  sub- 
domains  is: 

_  9T  ,  a2r  ,  a2T  .  . 

pCpdt  =  kx_Y0X2  +  kx~YQY2  +  Q  “  Grrans  (4) 


2.2.1.  The  current  collector  charge  balance  where  1  is  the  temperature  of  the  sub-domain,  p  is  the  density  of 

A  schematic  for  the  planar  and  transverse  current  density  the  electrode  plate  pair,  Cp  is  the  specific  heat  capacity  of  the 

distribution  on  the  two  current  collectors  including  the  tabs  is  electrode  plate  pair,  fex  -y  is  the  thermal  conductivity  of  the 


Tab  sub- domains 


Fig.  1.  The  computational  sub-domains  in  an  electrode  plate  pair. 
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Qj  =  -r-.(%j+ilYj)  (7) 

°CCJ  x  7 

If  the  transverse  heat  flux  follows  the  Newton’s  law  of  cooling, 
the  term  QTrans  is  expressed  as: 


Q-Trans 


h(T  —  Tamb) 
0.5  d 


(8) 


where  h  is  the  heat  transfer  coefficient,  d  is  the  distance  between 
the  two  planar  heat  transfer  surfaces,  and  Tam b  is  the  ambient 
temperature.  The  expressions  for  d  in  the  different  sub-domains  are 
listed  in  the  bottom  row  of  Table  2 


2.3.  The  models  in  the  mesoscopic  sub-domain 

The  electrochemical  processes  and  transport  phenomena  in  the 
porous  electrode  and  particles  are  described  by  the  pseudo-2D 
model  derived  by  Doyle  et  al.  [9]. 


2.3 A.  The  solution  phase  mass  transfer 

The  governing  equation  for  diffusion  of  lithium  ion  in  the 
solution  phase  is 


Insulation  boundary  conditions  are  assumed  at  the  interface 
between  the  coatings  and  the  current  collectors: 


Fig.  2.  Schematic  for  the  charge  balance  (current  flow)  in  current  collectors  and  tabs 
(discharge  process). 

electrode  plate  pair  in  the  X—Y  plane,  Q  is  the  volumetric  heat 
source,  and  Qirans  is  the  heat  transfer  rate  between  the  electrode 
plate  pair  and  the  surroundings. 

In  Equation  (4),  the  heat  source  term  Q  is  divided  into  three 
parts: 

Q  =®  £eQECh  +  £cc,pQ/,p  +  £cc,pQ/,n  (5) 

where  Qecii  is  the  electrochemical  heat  generated  in  the  porous 
electrode  sub-domain,  Qj,p  and  Q/  n  are  the  Joule  heating  rates  in  the 
current  collectors  and  tabs,  ze ,  £CCtP,  and  £Cc,n  are,  respectively,  the 
volume  ratio  factors  for  the  porous  electrode,  positive  electrode 
current  collector,  and  negative  electrode  current  collector.  The 
expressions  for  re,  £CCtP,  and  £Cc,n  in  different  sub-domains  are  listed 
in  Table  2.  Fig.  3  shows  the  thermal  conditions  for  a  unit  cross  area 
(in  X—Y  plane)  of  electrode  plate  pair.  The  planar  current  densities 
along  the  X—Y  plane  are  expressed  as: 


dce 


x=0 


=  0  and  — ^ 

ax 


=  o 


X — Ip+Is+In 


(10) 


where  aj  is  the  specific  surface  area  of  porous  electrode,  ce  is  the 
lithium  ion  concentration  in  the  solution  phase,  D2,eff  is  the  effec¬ 
tive  diffusion  coefficient  of  solution  phase,  t+  is  the  transference 
number,  Jj  is  the  rate  of  electrochemical  reaction  at  the  solid/solu¬ 
tion  surface,  and  ej  is  the  porosity  of  porous  electrode  j.  The  second 
term  on  the  right-hand  side  (the  source  term)  of  Equation  (9)  is 
zero  in  the  separator  since  there  is  no  reaction  in  that  region,  and 
the  porosity  £j  {j  =  p,n)  on  the  left-hand-side  is  then  replaced  by  the 
separator  porosity  es.  The  effective  diffusion  coefficient  in  solution 
phase  is  calculated  as: 

D2eff-=  DeMlkf>  j  =  p,s,n  (11) 


where  De,buik  is  the  diffusion  coefficient  in  bulk  solution  phase  (see 
Appendix  A),  and  (3j  is  the  Bruggeman  factor. 


a®/  .  .  efj  ,. 

1 XJ  =  and  iyj  =  -Occj-QY  0  =  P,n) 


(6) 


2.3.2.  The  solution  phase  charge  balance 

The  solution  phase  current  density  is  defined  as: 


The  Joule  heating  rates  in  the  Equation  (5)  are  expressed  as: 


Table  2 

The  expressions  for  the  volume  fraction  and  the  thickness. 


Parameter 

Sub-domain 

Tab  of  positive 
electrode 

Tab  of  negative 
electrode 

Cell 

0 

0 

lp  +  ls  +  In 

dp  +  Ip  +  h  +  In  +  dn 

£cc  ,p 

1 

0 

dp 

dp  +  Ip  +  Is  +  ln  +  dn 

£cc  ,n 

0 

1 

dn 

dp  +  Ip  +  ls  +  ln  +  dn 

d 

<5p 

<5n 

dp  +  Ip  +  Is  +  In  +  dn 

h 


2RT 


o 


alnce 

Ox 


(12) 


where  /ce ff  is  the  ionic  effective  electrical  conductivity  of  solution 
phase,  R  is  the  universal  gas  constant,  F  is  the  Faraday’s  constant, 
and  f±  is  the  activity  coefficient  of  solution  phase.  The  effective 
electrical  conductivity  of  solution  phase  is  calculated  as: 

zceff  =  Kbuik £jj  j  =  n  (13) 

where  /Cbuik  is  the  electrical  conductivity  in  bulk  solution  (see 
Appendix  A).  The  conservation  of  charge  in  the  solution  phase  is 
described  as 
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Fig.  3.  Schematic  for  the  thermal  conditions  of  cell. 


3  *2  CT 

to  =  a'FJ'  J  =  p  n 
with  boundary  conditions 

>2lx=o  =  0  and  *2  \x=in+is+ip  =  0 


(14) 


(15) 


h  lx=o  —  —  !'w.p 

(21) 

h\x=lp+ls+l„  =  {N,n 

(22) 

In  this  work,  the  one-direction  transverse  current  density  is 
defined  as: 


2.3.3.  The  solid  phase  charge  balance 

The  solid  phase  current  density  is  defined  as: 

h  =  J  =  P.n  (16) 

The  effective  electrical  conductivity  of  solid  phase  oy  eff  is  defined 


ft 

°fte ff  -  aj£sj 


(17) 


where  oj  is  the  electrical  conductivity  of  pure  active  material  for 
electrode  j,  and  eSj  is  the  volume  fraction  of  active  material  in  the 
electrode  j.  There  is  no  solid  phase  current  in  the  separator  region. 
The  conservation  of  charge  in  solid  phase  charge  is  described  as: 

§=  -ajFfj  j  =  P,n  (18) 


with  boundary  conditions: 

0i  l*=o  =  c^p-  h  \x=ip  =  h  \x=ip+is  = 

and  0i  lx=/p+/s+/n  =  (19) 


*N  =  *N,p  =  — i]V,n  (23) 

2.3.4.  The  electrochemical  reaction  kinetics 

The  surface  electrochemical  reactions  for  lithium  intercalation/ 
deintercalation  are  expressed  as  follow: 

Li  -  #s  <=±  Li+  +  e~  +  0S  (24) 


where  6S  represents  a  vacant  host  on  the  solid  particle  surface.  The 
rate  for  this  reaction  is  described  by  the  Butler-Volmer  equation 


Jj  —  ^rjce‘5(cs, 


\  0.5 
j,surf )  ( 


.0.5 
sj,  surf 


exp 


-  exp 


0.5F  \ 

~rTVj) 


(25) 


where  krj  is  the  reaction  rate  constant,  and  the  electrolyte 
concentration  in  solution  phase,  cSjiS urf  is  the  concentration  of 
lithium  at  the  surface  of  solid  phase  particles,  and  cSjimax  is  the 
maximum  lithium  concentration  in  the  solid  phase  particles.  The 
electrode  overpotential  r]j  is  defined  as: 


Vj  =  fa  -  <t>2  -  Uj  j  =  P,n 


(26) 


where  Uj  is  the  open  circuit  potential  which  depends  on  both  the 
surface  state-of-charge  and  temperature  (see  Appendix  A). 


Equation  (19)  suggests  that  the  solid  phase  charge  balance  in  the 
porous  electrode  sub-domain  is  coupled  to  the  macroscopic  cell 
sub-domain  at  the  two  boundaries  of  x  dimension  where  the 
potential  of  coatings  equals  to  the  potential  of  current  collector.  It 
can  be  derived  that  the  boundary  values  for  the  solid  phase  current 
density  at  the  two  current  collectors  are  equal: 

ill  x=lp+ls+ln  =  h\x=o  (20) 


2.3.5.  The  electrochemical  heat  source 

The  electrochemical  heat  source  term  QECh  in  Equation  (5)  is 
expressed  by  the  integrating  the  addition  of  volumetric  reversible 
heat  and  irreversible  heat  along  thickness  of  porous  electrode  sub- 
domain. 

Lp  \L>s  ~r L>n 

1 


QECh 


Ln 


I 


Jdx 


(27) 


The  terms  i^j  (where  j  =  p,n)  in  the  cell  sub-domain  charge 
balance  Equation  (1)  are  coupled  to  the  porous  electrode  sub- 
domain  through  the  following  correlations 


where  q'rev  and  q-rrev,  respectively,  denote  the  local  reversible  heat 
and  irreversible  heat  along  the  dimension  x,  and  are  defined  as 
follow: 
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<7rev  =  aJFJjT  ~Qf  J  =  P’n 

,  .  d(j>i  .  d(f>2  r, 

Cev  =  - 12 -g^  +  ajFJjVj  J  =  P,n 

As  the  overpotential  is  defined  as: 


(28) 


(29) 


QECh 


ip+is+in 

-iN  (4>p  -  4>„)  +  J  ajFJj  (fref-Qj  -  Uj, ref)  d* 

_ 0 _ 

Lp  +  Ls  +  ln 


(38) 


Let 


WE 


-iN(4>p-4>n) 

Fp  +  Ls  +  Ln 


(39) 


Vj  =  fa  ~  <t>2  -  Uj  j  =  p,n  (30) 

Equation  (29)  can  be  rewritten  as: 

Cev  =  ~  h  ^  +  ajFJjh  ~  ajFJj4>2  -  ajFJjUj  (31 ) 

Substitute  Equations  (14)  and  (18)  into  Equation  (31 )  and  obtain: 


and 

ip+is+in 

‘fc-VTC+C  /  <«» 

0 

And  the  heat  source  term  can  be  expressed  as  the  addition  of  WE 
and  Qh 

QEch  =  W£  +  Qw  (41) 


(32) 


Equation  (32)  can  be  rewritten  as: 

9(«101+«202) 


ax 


-Wj 


(33) 


Integrate  the  right-hand-side  of  Equation  (33)  from  0  to 
Ip  +  k  +  fn  along  the  x-dimension 


ip+/s+In  ip+4+in 

J  q;rrevdx=-(i^1+l202) J  cijFJjUjdx  (34) 
0  0 

According  to  the  boundary  conditions  of  charge  balance  in  the 
solution  phase  and  solid  phase  presented  in  Equations  (15)  and 
(19),  and  applying  the  correlations  shown  in  Equations  (20)— (23), 
Equation  (34)  can  be  simplified  to 


where  WE  describes  the  local  electrical  work  on  the  porous  elec¬ 
trode  and  Qh  describes  the  local  reaction  enthalpy  change 
throughout  the  porous  electrode  sub-domain.  According  to  the 
Thermal  Equations  (4)  and  (5),  the  energy  balance  in  the  cell  sub- 
domain  is  coupled  to  the  porous  electrode  sub-domain  through 
QECh-  For  the  high  rate  (i.e.  5C)  results  presented  in  this  paper,  we 
have  chosen  to  neglect  the  entropy  term,  dUj/dT ,  in  Equation  (40). 
This  term  should  be  included  for  low  rate  operations. 


2.3.6.  The  solid  phase  diffusion 

The  diffusion  of  lithium  in  the  solid  phase  follows  Fields  second 
law  in  a  spherical  coordinate  system,  and  the  governing  equation  is: 


5c, 


1  5 


at 


.2  ®Csj 


^  =  D^r¥i  y  =  p’n) 


(42) 


where  csj  is  the  concentration  of  lithium  in  the  solid  phase  parti¬ 
cles,  DSj  is  the  solid  phase  diffusion  coefficient  of  lithium.  The 
boundary  conditions  for  the  solid  phase  diffusion  are: 


-D  dCst 
usj  0r 


=  0 


r= 0 


(43) 


ip+is+in 

J  Qirrev^*  =  ~^n{^p  ~  ^n) 

0 


Ip+ls+ln 

J  cijFJjUjdx 

0 


(35) 


Substitute  Equations  (28)  and  (35)  into  (27),  so  that  the  heat 
source  term  can  be  expressed  as 


QECh 


lp+ls+ln 

-<N(*P  ‘M  ••  /  ajfJj  -  Uj) dx 

_ 0 _ 

Fp  +  Ls  +  Ln 


(36) 


The  temperature-dependent  open  circuit  potential  is  expressed 
in  Taylor’s  expansion: 


dUi 


Uj  =  Ujjef  [F  —  Tref 


(37) 


where  Tref  is  the  reference  temperature  and  UjtTef  is  the  open  circuit 
potential  at  reference  temperature.  Substitute  Equation  (37)  into 
Equation  (36)  and  obtain 


-D  •— s2 
SJ  dr 


r=Rs  i 


=  h 


(44) 


where  Rsj  is  the  radius  of  particle.  The  state  of  charge  for  electrode, 
6j,  is  defined  as: 

6j  =  -fU-  (45) 

csj,  max 

In  Equation  (26),  the  open  circuit  potential  of  electrode,  Uj,  is 
dependent  on  state  of  charge  6j  through  the  correlations  presented  in 
Appendix  A.  The  surface  lithium  concentration  cSj, SUrf  is  defined  as: 

csj, surf  =  csj\r=Rsj  (46) 

According  to  Equations  (25),  (44)— (46),  the  particle  sub-domain 
and  porous  electrode  sub-domain  are  coupled  through  cSjiSUrf  and  Jj. 


2.4.  The  temperature  dependency  of  parameters 

The  transport  and  kinetic  parameters  of  cell  materials  are  affected 
by  temperature.  The  diffusion  coefficient  and  ionic  electrical 


Table  3 

Parameter  values  for  LG  pouch  cell. 
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Sub-domain 

Parameter 

Value 

Positive  electrode 

Negative  electrode 

Particle 

Maximum  Li  capacity  csjjmax  (mol  m-3) 

49,000 

28,700 

Radius  of  particle  Rsj  (m) 

5.0  x  10-6 

12.5  x  10-6 

Reference  particle  diffusion  coefficient  Dsjref  (m2  s-1) 

3.0  x  10-15 

9.0  x  10-14 

Activation  energy  for  particle  diffusion  Eadij  (J  mol-1) 

4.0  x  103 

2.0  x  104 

Positive  electrode 

Separator 

Negative  electrode 

Porous  electrode 

Average  electrolyte  concentration  ce  (mol  m-3) 

1200 

Activity  coefficient  of  solution  phase /± 

0 

Thickness  lp  ls  ln  (m) 

50  x  10-6 

25  x  10-6 

70  x  10-6 

Bruggeman  factor  /fy 

2.0 

2.0 

2.0 

Porosity  rp  es  zn 

0.4 

0.4 

0.4 

Volume  fraction  of  active  material  eSj 

0.41 

0.51 

Specific  surface  area  ap  an  (m2  m-3) 

2.66  x  105 

1.02  x  105 

Solid  electronic  capacity  crp  an  (S  m-1) 

10 

100 

Reference  reaction  rate  constant  kr,p,ref  kr,n,ref  (m25  mol-0  5  s-1) 

4.966  x  10-11 

7.773  x  10-10 

Activation  energy  for  reaction  Eare,p  Eare,n  (J  mol  :) 

The  universal  gas  constant  R  (J  mol-1  K-1) 

3.0  x  104 

8.314 

3.0  x  104 

The  Faraday  constant  F  (C  mol-1) 

96,487 

Positive  electrode 

Negative  electrode 

Cell 

Current  collector  thickness  5j  (m) 

Current  collector  electronic  conductivity  accj  (S  m-1) 
Volumetric  heat  capacity  pCP  (J  K-1  m-3) 

Planar  thermal  conductivity  fcx-y(W  K  1  m-1) 

10  x  10-6 

37.8  x  106 

2.575  x  106 

27 

15  x  10-6 

59.6  x  106 

Positive  electrode 

Negative  electrode 

Tabs 

Current  collector  thickness  <5j  (m) 

Current  collector  electronic  conductivity  accj  (S  m-1) 
Volumetric  heat  capacity  pCP  (J  K  1  m-3) 

Planar  thermal  conductivity  kx~ y  (W  K  1  m-1) 

10  x  10  6 

37.8  x  106 

2.42  x  106 

237 

15  x  10-6 

59.6  x  106 

3.41  x  106 

401 

conductivity  of  solution  phase  are  dependent  on  temperature  and 
electrolyte  concentration  through  the  empirical  correlations  pre¬ 
sented  in  Appendix  A.  The  diffusion  coefficients  of  solid  phase  and  the 
rate  constants  for  the  intercalate/deintercalate  reactions  are  depen¬ 
dent  on  temperature  through  Arrhenius’  equations: 


Dsj  —  Dsj  ref£Xp 


Eadij  ! 

|T_ 

lV 

R  1 

V 

Tref  J 

(47) 


krj,refexp 


Earej  _J_A 

R  [T  Tre[J 


(48) 


where  DSJiref  and  /<rjiref  are  respectively  the  solid  phase  diffusion 
coefficients  and  reaction  rate  constants  at  reference  temperature  Tre f, 
and  Eadij  and  Earej  are,  respectively,  the  activation  energy  for  solid 
phase  diffusion  and  intercalate/deintercalate  reactions  (see  Table  3). 


Input  variables : 


Output  variables : 


Node  point  1  in  Ar-dim  Node  point  i  in  Ar-dim 


Fig.  4.  The  pseudo-2D  model  expressed  as  control  block. 
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Node  point  1  in  the 
cell  sub-domain 


Node  point  2  in  the 
cell  sub-domain 


Node  point  i  in  the 
cell  sub-domain 


a 


Cell  and  Tab  sill 

_  &F  ,  d2T  ,&T 

s2®  a2o„  if 

c  p  +  cr  p  +  N 

ttJf  3Z2  3F2  0.55, 

f{t)~llTdXdr  *.(t) 

b-domains  (Z  and  Y  di 

2jj>  +  S,QeCW  +  Scc  -  £ 

( <J>.  -  )  +  a^Brh 

_ V  ?  >  94> 

r  * 

A  c 

-  qj^2  *  °ic, 

\dXdY  « 

A:,]1 

imensions) 

?Trms 

(*„-*») 

<*>p, 

*  =0 

ar2  0.54  4:.i 

;do~ii 

^Cell 

&*(«*.  r) 

_  _ 
J 

ELectro  chcnii  c  al 

T(i\ 

model  block 

b 


Fig.  5.  The  two  types  of  model-coupling  approaches:  (a)  the  full-distribution  model,  (b)  the  linear  model. 


2.5.  The  electrical  limiting  condition 

According  to  Equation  (3),  there  is  no  imposed  boundary  condition 
for  the  current  collector  charge  balance;  therefore,  two  additional 
equations  are  needed  to  determine  the  unique  values  for  the  electrical 
potentials  of  the  current  collectors.  The  first  equation  is  obtained  by 
setting  the  potential  at  the  top-center  of  negative  electrode  tab  to  be 
zero,  and  the  second  equation  is  the  electrical  limiting  equation  that 
the  integration  of  transverse  current  density  through  the  cell  sub- 
domain  equals  the  total  applied  current  on  the  electrode  pair: 

J  j  iNdXdY  =  -/app  (49) 

^Cell 

where  Japp  is  the  total  applied  current  passing  through  the  electrode 
pair  and  is  defined  as  negative  for  discharge. 


2.6.  The  coupling  of  models  of  different  length  scales 

In  this  work,  the  Newman’s  pseudo-2D  model  in  the  mesoscopic 
dimensions  was  expressed  as  a  control  block  as  presented  in  Fig.  4. 
The  inputs  for  the  control  block  are  the  temperature  T  which  affects 
certain  transport  and  kinetic  processes,  and  the  current  collector 
potentials  4>j  ( j  =  p,n)  which  act  as  boundary  values  for  the  solid- 
phase  charge  balance  at  the  interfaces  between  coatings  and 
current  collectors.  The  outputs  from  the  control  block  are  the 
transverse  current  density  z'n  and  the  electrochemical  heat  source 
QECh-  As  the  input  variables  are  distributed  with  the  macroscopic 
dimensions,  the  pseudo-2D  control  block  has  to  be  solved  at  every 
node  point  of  the  cell  sub-domain  and  then  return  the  local  elec¬ 
trical  and  thermal  outputs  (Fig.  5(a))  to  the  charge  balance  Equation 
(1)  and  the  energy  balance  Equation  (4).  In  this  work,  this  model 
coupling  approach  as  described  in  Fig.  5  is  named  as  the  full- 
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45.0 


45.0 
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145.0 


(mm) 


Fig.  6.  The  dimensions  of  electrode  in  the  LG  pouch  cell. 


distribution  model.  The  disadvantage  of  full-distribution  model  is 
that  a  very  huge  non-linear  differential-algebraic  equation  (DAE) 
system  is  to  be  solved  in  simulations.  Therefore,  the  full- 
distribution  model  is  very  time-consuming  although  it  accurately 
describes  the  non-uniformity  of  cell’s  thermal  and  electrical 
behaviors. 

Alternatively,  in  our  work,  a  simplification  method  is  proposed 
to  produce  our  ROM  and  reduce  the  simulation  time  and  memory 
requirements.  Instead  of  using  the  reduced-order  model  to  replace 
pseudo-2D  model,  our  simplification  method  focuses  on  the 
segregation  of  sub-domains.  As  shown  in  (Fig.  5b),  the  distributed 
input  variables  for  the  pseudo-2D  control  block  are  replaced  by  the 
volume-average  of  those  variables  through  the  cell  sub-domain; 
the  electrochemical  heat  source  returned  to  the  macroscopic  sub- 
domains  are  approximated  linearly  through  the  following  Taylor’s 
expansion: 


QECh~QECh(^p5  ^n, 


9QECh 


d<f>p 

9QECh 


d$n 


*n,  T 

T 


(®p  -  ®p) 

($n  -  ®n) 


(50) 


Equation  (50)  was  derived  by  assuming  that  the  spatial  varia¬ 
tion  of  current  collector  potential  is  the  major  reason  for  the 
distribution  of  the  electrochemical  heat  source.  This  assumption 
can  be  verified  by  Equation  (38)— (41),  Equation  (39)  suggests  that 
the  electrical  work  on  the  porous  electrode  is  almost  a  linear 
function  of  the  current  collector  potentials,  where  the  derivatives 
of  We  to  4>j  are: 


dWg  = _ jjy  anc|  dW£  =  iN 

d4>p  Lp  +  Ls  +  Ln  Lp  +  Ls  +  Ln 


Case  1 


Case  2 


Newton  cooling 

h  =  lW-m^-K-1 

for  all  boundary  and 
transversal  heat  flux 


Case  3 


Fig.  7.  The  boundary  and  transversal  thermal  condition  settings  for  different  cases. 
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Table  4 

The  computation  time  for  case  studies. 


Model  type 

Number  of  equations 

Simulation  time 

Case  1 

Case  2 

Case  3 

Full-distribution 

7576 

5510  s 

6371  s 

5498  s 

Linear 

361 

6.26  s 

6.33  s 

6.60  s 

Equation  (40)  suggests  that  the  reaction  enthalpy  Qh  only 
depends  on  the  entropy  change  and  the  electrode  open  circuit 
potential  at  the  reference  temperature,  and  the  distributions  of 
these  variables  in  the  macroscopic  dimensions  are  neglected.  The 
transverse  current  returned  from  the  pseudo-2D  control  block  was 
replaced  by  the  volume-average  value  in  the  cell  sub-domain: 


a 


>- 


Temperature  (°C) 


_ XJm] _ 

| Full-distribution  model 


Temperature  (°C) 


Fig.  8.  The  simulated  end-of-discharge  temperature  profiles  for  different  cases:  (a) 
case  1,  (b)  case  2,  (c)  case  3. 


In-  ~ 


Japp 

Acell 


(52) 


where  Aieii  is  the  area  of  the  cell  sub-domain.  In  this  work,  this 
approximated  model  coupling  approach  is  called  as  the  linear  model; 
and  through  this  method,  the  pseudo-2D  electrochemical  model  was 
only  evaluated  once  at  the  average  temperature  and  potential  in  the 
cell  sub-domain,  time  and  computer  memory  was  significantly  saved. 


3.  Results  and  discussion 

In  proposing  the  linear  approximation  approach,  our  concerns 
include  how  much  computation  time  can  be  saved  by  decoupling 
the  sub-domains  and  how  well  the  linear  model  can  agree  with  the 
full-distribution  model.  Case  studies  were  also  made  in  this  work  to 
analyze  the  distribution  of  heat  source  and  planar/transverse 
current  density.  As  our  model  was  just  based  on  one  single  planar 
electrode  pair,  the  model-to-data  validations  will  be  made  in  future 
work  until  the  model  is  extended  to  the  cell  stacks  with  multiple 
electrode  pairs. 

In  our  work,  the  model  was  based  on  the  electrode  plate  pair  of 
LG  pouch  cell  (see  Kim  et  al.  [7]),  the  electrode  dimensions  are 
presented  in  Fig.  6.  The  physical  and  chemical  properties  of  the 
electrodes  are  presented  in  Table  3,  where  the  active  materials  for 
the  positive  and  negative  electrodes  are,  respectively,  LiNi0.8- 
C00.15AI0.05O2  (NCA)  andLi*C6.  This  model  was  solved  using  the 
Galerkin’s  finite  element  method;  the  setting  of  geometry,  the 
meshing  of  sub-domains,  the  solving  of  DAE  system,  and  the  post¬ 
processing  of  results,  were  all  performed  by  MATLAB. 

The  comparisons  between  the  linear  model  and  the  full- 
distribution  model  were  made  by  simulating  the  5C  (where 
C  =  0.4186  A  for  one  electrode  pair)  discharge  of  electrode  with 
three  different  boundary  and  transversal  thermal  conditions  as 
presented  in  Fig.  7.  The  initial  state-of-charge  (SOC)  is  0.41  for 
positive  electrode  and  0.63  for  negative  electrode,  the  cut-off 
voltage  is  3  V,  and  both  the  ambient  temperature  and  the  initial 
temperature  of  cell  are  25  °C.  The  computation  time  that  the  linear 
model  and  the  full-distribution  model  take  for  each  simulation  are 
listed  in  Table  4,  and  the  linear  model  is  faster  than  the  full- 
distribution  model  by  about  900  times. 

The  temperature  profiles  at  the  end  of  5C  discharge  simulated  by 
the  linear  and  the  full-distribution  models  for  the  three  different 
cases  are  presented  in  Fig.  8(a)— (c).  According  to  these  plots,  the 
maximum  temperature  difference  between  the  hottest  and  the 
coldest  spots  of  the  cell  is  16.5  °C  (in  case  2);  and  in  these  cases,  the 
linear  model  shows  excellent  agreements  with  the  full-distribution 
model.  The  results  in  Fig.  8  also  show  that  the  temperature  of  the 
positive  electrode  tab  is  always  higher  that  that  of  the  negative 
electrode  tab,  the  reason  is  that  the  current  collector  of  positive 
electrode  which  is  made  by  the  aluminum  foil  has  larger  electrical 
resistivity  than  the  current  collector  of  negative  electrode  which  is 
made  by  the  copper  foil,  and  thus  more  Joule  heat  is  generated  on 
the  positive  electrode  than  on  the  negative  electrode.  In  Fig.  8(b) 
and  (c)  (case  2  and  case  3),  the  temperature  of  the  two  electrode 
tabs  is  lower  than  that  of  the  cell  sandwich  when  transversal  and 
boundary  thermal  conditions  are  same  for  both  cell  and  tabs.  The 
reason  is  that  the  tabs  have  much  smaller  thickness  than  the  cell 
sandwich;  and  according  to  Equation  (8),  the  transversal  heat 
transfer  rate  is  inversely  proportional  to  the  thickness,  therefore 
the  heat  dissipates  much  faster  on  the  tabs  than  on  the  surface  of 
cell  sandwich.  In  Fig.  8(a)  (case  1)  where  the  heat  transfer  coeffi¬ 
cient  on  the  tabs  is  smaller  than  that  on  the  cell  sandwich,  the 
negative  electrode  tab  is  the  coldest  part;  the  hottest  spot  locates  in 
the  top-left  part  of  the  cell  sub-domain  which  is  close  to  the  tab  of 
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Point  2 


Point  1 

J 


Point  3 


Fig.  9.  The  temperature  vs  time  plots  at  different  points  with  thermal  condition  case  1. 


the  positive  electrode;  and  temperature  of  other  locations  in  the 
cell  sun-domain  lies  between  the  temperature  of  the  positive  and 
negative  electrode  tabs.  The  plots  for  variation  of  temperature  at 
three  selected  points  with  time  during  case  1  are  presented  in 
Fig.  9,  and  the  linear  model  still  agrees  well  with  the  full- 
distribution  model.  As  shown  in  Fig.  9,  the  temperature  of  the 
upper  point  (point)  rises  rapidly  at  the  beginning  (0-100  s)  and 
then  the  increase  of  temperature  slows  down;  the  temperature  at 
lower  points  (points  1  and  2)  just  increases  with  nearly  constant 
rates.  The  reason  is  that  the  temperature  rise  at  point  3  is 
decelerated  by  the  high  thermal  dissipation  rate  on  the  tab  of 
positive  electrode. 


Fig.  10.  The  end-of-discharge  potential  profiles  with  thermal  condition  setting  case  1. 


The  case  studies  above  also  show  the  great  potential  values  for 
this  model; 

i)  This  model  could  be  extended  for  life  study  by  including 
capacity  fade.  The  results  suggest  that  the  high  non¬ 
uniformity  of  the  cell  thermal  behavior  could  significantly 
affect  the  cycling  life  of  Li-ion  cells. 

ii)  This  model  could  also  be  used  for  the  design  and  optimization 
of  cells  and  packs.  For  example,  the  results  above  show  that 
the  tab  temperature  is  strongly  affected  by  the  thickness,  and 
the  electrical  properties  of  the  tab  material.  This  model  could 
be  used  to  find  appropriate  tab  dimensions  to  minimize  the 
high  rate  tab  temperatures. 

The  current  collector  potential  profiles  at  the  end  of  5 C 
discharge  simulated  by  the  linear  and  the  full-distribution  models 
for  thermal  condition  case  1  are  presented  in  Fig.  10(a)  and  (b),  and 
the  linear  model  still  agrees  very  well  with  the  full-distribution 
model.  As  shown  in  Fig.  10,  the  electrical  potentials  of  current 
collectors  change  by  10-20  mV  through  the  cell  and  tab  sub- 
domains,  and  the  potential  gradients  in  the  tabs  are  larger  than 
that  in  the  cell  sub-domain  due  to  the  high  planar  current  density 
in  the  tabs. 


4.  Conclusion 

In  this  work,  a  multi-scale  multi-dimensional  thermal  model  for 
Li-ion  cell  was  developed  and  a  linear  approximation  was  also 
proposed  to  simplify  this  model.  The  linear  model  is  much  more 
efficient  than  the  full-distribution  model  in  terms  of  the  compu¬ 
tation  time.  The  linear  model  also  maintains  high  accuracy  in  the 
predictions  of  most  model  variables  except  the  two  coupling  vari¬ 
ables  which  are  directly  approximated  linearly  to  decouple  the  sub- 
domains,  but  this  error  does  not  affect  the  evaluation  of  other 
variables.  Therefore,  the  linear  model  is  an  effective  approach  to 
study  the  distributed  electrical  and  thermal  behavior  of  an  elec¬ 
trode  plate  pair  in  a  Li-ion  cell.  In  the  future,  this  model  will  be 
extended  to  include  multiple  electrode  plate  pairs  and  model 
predictions  will  be  further  validated  by  experimental  data. 


Appendix  A 

The  physical  properties  of  the  solution  phase  are  shown  in 
Equations  (Al)— (A3).  In  these  equations,  the  solution  phase 
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concentration  ce  is  in  mol  m  3  and  temperature  T  is  in  I<.  The 
diffusion  coefficient  of  bulk  electrolyte,  De,buiio  is  of  the  unit  m2  s_1. 


D, 


e,bulk 


=  5.84  x  10  7exp(  - 


/  2870\ 

(  Ce  ) 

V  T  ) 

VI 0007 

-33.9 


xl0-7exp(-^°)(Ig5)+129 


x 10  7exp(  - 


3200\ 

T  J 


(Al) 


The  electrical  conductivity  of  bulk  electrolyte,  Kbuik.  is  of  the  unit 


S  m 


e“k  - 

+244axp(-^2)(T^5)  (A2) 

The  transference  number  of  electrolyte,  t+,  is  dimensionless. 


t+  =  2.67  x  10-4exp(^)  (t^q)2+3  09 


x  10  3exp 


/653\ 

(  Ce  ) 

{  T  ) 

viooo; 

i.517exp 


(A3) 


The  open  circuit  potentials  of  the  negative  and  positive  elec¬ 
trodes  are  shown  in  Equations  (A4)  and  (A5),  where  Un  and  Up 
are  in  V. 


U„  =  0.124  +  1.5exp(— 7O0„) 

-0.0351  tanh(12.0482  dn  -  3.4458) 

-  0.0045  tanh(8.4034  8„  -  7.5630) 

-  0.035  tanh(200  8n  -  19.80) 

-  0.0147  tanh(29.4118  6n  -  14.7059) 

-  0.102  tanh(7.0423  6„  -  1.3662) 

-  0.022  tanh(60.9756  dn  -  59.7561) 

-  0.01 1  tanh(44.2478  d„  -  5.4867) 

-0.0155  tanh(34.4828  dn  -  3.6207)  (A4) 


Up  =  1.638  8™  -  2.222  d9p  +  15.056  d8p  -  23.488  67p 

+  81 .246  8 p  -  344.566  8Sp  +  621 .3475  84p  -  554.774  8p 
+  264.427  82p  -  66.3691  0P  + 11 .8058 
-0.61386  exp  (5.8201  6 »’36'4) 

(A5) 
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List  of  symbols 

ay  Specific  surface  area  of  electrode  (m2  m-3) 

Are//-'  Area  of  the  cell  Sub-domain  (m2) 

ce:  Concentration  of  Li  in  solution  phase  (mol  m-3) 

ce  :  Average  electrolyte  concentration  (mol  m-3) 

Cp:  Specific  heat  capacity  (J  K-1  kg-1) 

csy  Concentration  of  lithium  in  the  solid  phase  particles  (mol  m-3) 

Csj.max-  Maximum  Li  capacity  in  solid  phase  (mol  m-3) 
d:  Distance  between  the  two  planar  heat  transfer  surfaces  (m) 

D2,eff  Effective  diffusion  coefficient  of  solution  phase  (m2  s-1) 

De,buik'-  Diffusion  coefficient  in  bulk  solution  phase  (m2  s_1) 

Dsf  The  temperature-dependent  particle  diffusion  coefficient  (m2  s  1) 

Dsjre/-'  Reference  particle  diffusion  coefficient  (m2  s-1) 

Eadif.  Activation  energy  for  particle  diffusion  (J  mol-1) 

Earey  Activation  energy  for  reaction  (J  mol-1) 
f±:  Activity  coefficient  of  solution  phase 
F:  The  Faraday  constant  (C  mol-1) 

h:  Heat  transfer  coefficient  for  Newton’s  law  of  cooling  (W  K-1  m-2) 
if.  Solid  phase  current  density  (A  m-2) 
if.  Solution  phase  current  density  (A  m-2) 

iappf  Applied  outward  current  density  on  the  top  of  electrode  tabs  (A  m-2) 

Current  density  through  the  coating/current-collector  intersection  (A  m-2) 
if.  Transverse  current  density  (A  m-2) 
ixy  Planar  current  density  in  X  direction  (A  m-2) 
ixy  Planar  current  density  in  Y  direction  (A  m-2) 
lapp.  Total  applied  current  on  the  electrode  pair  (A) 

Jy  Rate  of  electrochemical  reaction  at  the  solid/solution  surface  (mol  m-2  s_1) 
krj:  The  Temperature-dependent  reaction  rate  constant  (m2  5  mol-0  5  s-1) 
krj.ref-  Reference  reaction  rate  constant  (m2  5  mol-0,5  s-1) 
kx-Y :  Planar  thermal  conductivity  (W  K-1  m-1) 

If  Thickness  of  electrode  j  or  separator  (m) 
n:  Unit  normal  vector  pointing  out  of  the  boundary 
q'rev:  Local  reversible  heat  (W  m-3) 
q'irrev:  Local  irreversible  heat  (W  m-3) 

Q-'  The  total  heat  source  (W  m-3) 

Qec/i-  The  electrochemical  heat  (W  m-3) 

Oh-  Reaction  enthalpy  change  throughout  the  porous  electrode  sub-domain 
(W  itT3) 

Qjf.  The  Joule  heat  in  the  current  collectors  (W  m-3) 

Qirans'-  Transversal  heat  transfer  rate  (W  m-3) 

R:  The  universal  gas  constant  (J  mol-1  K-1) 

Rs:  Radius  of  active  material  particles  (m) 
t+:  Transference  number  of  electrolyte 
T:  Temperature  (K) 

Tamb:  Ambient  temperature  (K) 

Tref ■  Reference  temperature  (K) 

Uj:  Open  circuit  potential  of  electrode  (V) 

We:  Electrical  work  on  the  porous  electrode  sub-domain  (W  m-3) 
ft:  Bruggeman  factor  of  electrode  or  separator 
ft:  Thickness  of  current  collector  j  (m) 

Ee:  Volume  ratio  factors  of  the  porous  electrode 
£CCJ-:  Volume  ratio  factors  of  the  current  collector 
Eji  Porosity  of  electrode  or  separator 

£SJ-:  Volume  fraction  of  active  material  in  the  porous  electrode 
7]f.  The  electrode  overpotential  (V) 
ft:  Electrode  state  of  charge 

Keff:  Ionic  effective  electrical  conductivity  of  solution  phase  (S  m-1) 

k bulk:  Ionic  electrical  conductivity  in  bulk  solution  (S  m-1) 

p:  Density  of  electrode  pair  (kg  m-3) 

ay.  Electrical  conductivity  of  pure  active  material  (S  m_1) 

ajiefy  Effective  electrical  conductivity  of  solid  phase  (S  m_1) 

0j:  Electrical  potential  of  solid  phase  (V) 

02-  Electrical  potential  of  solution  phase  (V) 

Oy  Electrical  potential  of  current  collector  (V) 

Qceii-  Symbol  of  the  cell  sub-domain  in  integral 


